Absence of gapped broken inversion symmetry phase of electrons in bilayer graphene 

under renormalized ring-diagram approximation 
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On a lattice model, we study the possible existence of a gapped broken inversion symmetry phase 
(GBISP) of electrons with long-range Coulomb interactions in bilayer graphene using both the self- 
consistent Hartree-Fock approximation (SCHFA) and the renormalized ring-diagram approximation 
(RRDA). The RRDA takes into account the charge density fluctuations beyond the SCHFA. Al- 
though the GBISP at low temperature and low carrier concentration is predicted by the SCHFA, 
we show here that this state can be substantially suppressed by the charge density fluctuations in 
the RRDA. We also present a numerical algorithm for calculating the self-energy of electrons with 
the singular long-range Coulomb interaction on the lattice model. 

PACS numbers: 71.10.-w,71.10.Fd,73.22.Pr,71.15.Dx 



I. INTRODUCTION 

Because of its tunable band gap, which can be changed 
through an external gate voltage, bilayer graphene is a 
promising material with a great potential for application 
to new electronic devicesii"— In the low-carrier-doping 
regime of bilayer graphene, electrons are strongly coupled 
via Coulomb interactions. The phase of bilayer graphene 
in this low carrier concentration, and at low tempera- 
ture, is still not completely understood. Several candi- 
dates have been suggested for the ground state, such as a 
ferroelectric-layer asymmetric state^^^— a layer-polarized 
antiferromagnetic state /"'"'^^ a quantum anomalous Hall 
state^ '^^1-^^ a quantum spin Hall state^ii^ a quantum 
valley Hall statej^ a charge density wave state^^ and 
the possibility of gapless states, such as the nematic 
statej ^^i^'' The experimental observations on the ground 
state of bilayer graphene, all performed on high quality 
suspended samples, are also controversial. Some experi- 
mental results showed that the system is gapped at the 
neutrality point, ^^^^^ whereas one experiment found a 
gapless statei^^ So far, most of the theoretical studies 
are based on the self-consistent Hartree-Fock approxi- 
mation (SCHFA), ^'^'^^ many-body perturbation theory,^ 
and the renormalization group approachJii^iii All the 
above approaches have been applied to the simplified two- 
^5,7,9 ,11,14,17 g^jj^ four-band^^ continuum models. It is well 
known that the SCHFA usually overestimates the order 
parameter characterizing a broken symmetry phase and 
the transition temperature because it neglects the fluc- 
tuations of the effective one-body interaction field and of 
other one-body observables such as the charge density. 
Since the understanding of the electronic state of bilayer 
graphene at low carrier doping and low temperature is a 
fundamental issue for graphene physics, it is necessary to 
investigate the state with a more sophisticated approach 
that takes into account the effect of charge density fluc- 
tuations on top of the mean-field ground state. 

In this work, we study the existence of a gapped bro- 
ken inversion symmetry phase (GBISP) using both the 



SCHFA and the renormalized ring-diagram approxima- 
tion (RRDA).^-^ The RRDA takes into account the charge 
density fiuctuation (CDF) effect beyond the mean field 
and satisfies the microscopic conservation laws.^"* For an 
electron system with long-range Coulomb interactions, 
CDF is the predominant contribution to the self-energy 
of electrons. It has been shown^^ that the RRDA results 
for the ground-state energy of two- and three-dimensional 
interacting electron gases are more accurate than the 
random-phase approximation (RPA) results when com- 
pared with Monte Carlo simulations. 

In the RRDA, the Green's function and self-energy 
are self-consistently determined by coupled integral equa- 
tions. The self-consistent calculation of the self-energy 
in momentum space involves carrying out many convo- 
lutions, which are numerically expensive. In order to re- 
duce the computational time required by our approach, 
we convert the convolutions in momentum space to mul- 
tiplications in real space. Since the continuum model 
is the low-energy limiting case of the lattice model, the 
momentum of the electrons is confined within two val- 
leys around the Dirac points i^^'^'' Because of the finite 
momentum cutoff for each valley, the conversion of the 
convolution from momentum space to real space is no 
longer valid for the two- and four-band continuum mod- 
els. Instead of modeling bilayer graphene with an effec- 
tive continuum model, we therefore sketch it as a bilayer 
of a hexagonal lattice model. The lattice model does not 
require a momentum-space cutoff, and is therefore im- 
mune to the aforementioned problems of the continuum 
models. 

The key problem in calculating the self-energy is to 
manage to deal with the long-range Coulomb interaction 
between electrons accurately. For the two-dimensional 
system under consideration, this interaction is inversely 
proportional to the momentum transfer q in the long- 
wavelength limit. In a continuum model, one can trans- 
form the 1 /q singularity to the logarithmic one after per- 
forming the azimuthal integration^^ and then get rid of 
the logarithmic singularity by special treatment. In a lat- 



tice model, however, we cannot perform the azimuthal in- 
tegration analytically and must face the 1/q singularity. 
Since dealing with the long-range Coulomb interaction 
is inevitable in many-body problems, we now present a 
numerical algorithm to tackle the interaction divergence 
issues systematically. 

II. LATTICE MODEL 

The lattice structure of bilayer graphene is shown in 
Fig. 1. The unit cell in each layer is represented by a 
diamond. The unit cell of the bilayer system contains 
four atoms denoted as ai, bi, a.2 and b2. The lattice 
constant of monolayer graphene is defined as the distance 
between two nearest corner atoms in the diamond and is 
given by a « 2.4 A . The interlayer distance is zq — 3.34 
A « 1.4a. The energy of electron hopping between the 
nearest-neighbor (NN) carbon atoms in each layer is t 
2.82 eV^ while the interlayer NN hopping is ti w 0.39 
eV^ 

The Hamiltonian describing the electrons is given by 
= - ^ Uj cL Cja + ]^Y1 '^"J ( 1 ) 

ija ij 

where creates an electron at site i with spin cr, 
Srij — rij — n with rij as the electron density operator 
at site J and n the average occupation number of elec- 
trons per site (which is also the charge number of the 
neutralizing background), and Vij is the Coulomb inter- 
action between electrons at sites i and j. The model 
is restricted to NN hopping within the same layer and 
between the adjacent sites on top and bottom layers as 
shown in Fig. 1. As described by Eq. ([T]), we con- 
sider here only the charge-charge interactions. Since the 
long-range antiferromagnetic order is prohibited'^- in two- 
dimensional space, we neglect the antiferromagnetic cou- 
pling due to the on-site repulsion in the present work. 

We now consider the behavior of Coulomb interaction 
Vij between two electrons at sites i and j. At long dis- 
tance, Vij is given by Vij = /erij with e the dielectric 
constant in the high-frequency limit of the system and r.^ 
the distance. However, at short distances, because of the 
spread of the 7r-orbital wave function of the conduction 
electrons, Vij is weakened from the behavior 1/rij. Tak- 
ing the effect of the wave function spread into account, 
we model the interaction as 

Vij = [1 - exp(-ry/ro)] (2) 

srij 

with ro = a. Clearly, Vij behaves as e'^/erij at large r^ , 
while it is suppressed from the 'bare' Coulomb interaction 
{e^ /cTij) at small rij. In particular, at — 0, it is given 
by a finite value e^/erp. For the present electron sys- 
tem with long-range Coulomb interactions, the final re- 
sult under consideration should not be sensitively depen- 
dent upon the details of the short-range behavior of the 




FIG. 1: Left: Structure of Bernal stacking bilayer graphene. 
Right: Top view of the bilayer graphene. The parameters t 
and ti axe the electron hopping energies between one atom 
and its nearest neighbor belonging to the same layer, and to 
the neighboring layer above or below, respectively. The unit 
cell of each layer is represented by the green-sided diamond. 



interaction. This can be understood from the behavior of 
its Fourier component in momentum space. The Fourier 
component is singular at the long-wavelength limit and 
the singular part is independent of the short-range be- 
havior. At low carrier concentration, the electrons state 
is mainly determined by the singular part of the inter- 
action. We use the dimensionless constant g = /eat 
to denote the strength of Coulomb coupling. The range 
0.4 < g < 1.8 covers the cases of various experimental 
setups, from suspended bilayer graphene (BLG) to BLG 
placed on substrates^ as Si02 and ice. 

The system defined by Eq. ^ satisfies the particle- 
hole symmetry. To see this, we denote the doped electron 
concentration per carbon atom as i5 and have n = \ + 5. 
Under the transformation 5 — > —5 and Cj^a — > +(— )cj ^ 
and c j ^ +(— )cj^cr for electrons at a,j [hj) sites, H is 

unchanged. Furthermore, K — H — ii{N — A^o) {N being 
the total electron number operator and Nq being the total 
number of lattice sites, so that the operator N — Nq refers 
to the total number of doped electrons) is also unchanged 
under the above electron-hole transformation, provided 
fi — !■ —fi. Thus, the chemical potential /i must be an odd 
function of S. 

The Green's function G of the electron system is de- 
fined as 

G(z, j,r - r') = -{TrQAr)Clir')) (3) 

where C}^ = (4^,,, 4^^-,, 4^-,, 4^-^) with c)^,(fc^)j-, cre- 
ating an electron of spin a at site a/ (b/) of the ^th (= 
1,2, respectively, for top and bottom) layer of the jth unit 
cell. In momentum- frequency space, G (a 4x4 matrix) 
can be expressed in terms of the self-energy T.(k, iui) as 

G{k, iuJi) = [iuji -t- - /ife - S(A:, iuie)Y^ (4) 
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FIG. 2: Self-energy of the SCHFA. Left: Hartree term. Right: 
Fock exchange term. The sohd hne with an arrow denotes the 
Green's function. The wavy hne is the Coulomb interaction. 



with 



hh = 




(5) 



where lu^ = (21 + 1)tiT is the fermionic Matsubara fre- 
quency with i as integer number and T the temperature, 
and e/c = —t[l + ex.p{—ikx) + exp(— ifc^)]. Here ^ is the 
chemical potential and is determined by 



2T 



— > TrG'(fc, iujt) exp{iujeri), 
l\o — 



(6) 



where the factor 2 stems from the spin degeneracy and rj 
is an infinitesimally small positive constant. To proceed, 
we need to provide an approximation for Yi{k,iuje). In 
the following sections, we investigate the possibility of 
the existence of the GBISP using the SCHFA and the 
RRDA for the self-energy, respectively. 



III. STUDYING THE EXISTENCE OF THE 
GBISP USING THE SCHFA 



Hartree term is diagonal, T,f[^ = IS. ^5^^, with 



A„ = 6iu 



pi 



Uf_,i = \im[Vf_,i{q) - Vf_,4{q)] 



-Up2 = lim[wp2('7) 



"^3 



(7) 
(8) 

(9) 



where Vf_iu{q) (with the subscripts being the same 
as those used in the definition of the Green's function, 
denoting the four sublattices ai, bi, a.2 and b2) is the 
Fourier component of the Coulomb interaction. In the 
long- wavelength limit, v^^{q) behaves like 



2^ 



exp(-z^j,(5) + v^,i,, q-)-0 (10) 



where 5*0 = \/3a'^/2 is the area of the two-dimensional 
unit cell of monolayer graphene, and Q is the magnitude 
of the vector Q — Mq with^^ 



M 



1 





1 2 

V3 V3 



(11) 



and where the components of q are along the nonorthog- 
onal axes of the diamond-shaped Brillouin zone. The 



value of 



or zq (the distance of the two layers) 



depends on fiv denoting the same layer or two different 
layers. The last term in Eq. (1101) . v^^, is the regular part 
of the Coulomb potential for g — > 0. The q dependence in 
Eq. (llOp is different from the conventional form because 
the coordinate axes of the reciprocal lattice where q is 
defined are nonorthogonal. The wave vector Q is defined 
in an orthogonal basis The relations Ai = — A4 and 
A2 = — A3 can be easily checked. 
The Fock exchange term is given by 



(12) 



Let us first consider the physical meaning of the 
GBISP. As can be seen in Fig. 1, supposing the ori- 
gin is at the middle point of a bia2 bond, when chang- 
ing each atom at site rj to —rj, the whole lattice is 
unchanged. This transformation is equivalent to inter- 
changing the top and bottom layers and then rotating 
the lattice by an angle tt around the bia2 bond. In the 
non symmetry-broken state, the electron system is un- 
changed with respect to such a transformation. However, 
when the strong Coulomb interactions drive the system 
to a GBISP, the two layers cease to be equivalent by in- 
version symmetry, and the electrons experience different 
fields on the two layers. Specifically, there may exist net 
electronic charge accumulation at each atom. We denote 
the deviations of the electronic charge density from the 
average value n at each of the four sites of the unit cell 
as {61,62,-62,-61). 

Under the SCHFA or the mean-field approximation, 
the self-energy is diagrammatically given by Fig. 2. The 



where M = iVo/4 is the total number of unit cells in one 
layer, and n^y{k) is given as 

npi/(fc) = G ^u[k,iui) e-yiY>{iujtTii) - 5^,^/2, (13) 



which corresponds to the quasiparticle distribution func- 
tion, the term —5^i^/2 stemming from the non-normal 
order of the electronic interaction operator. Under 
the mean-field approximation, the self-energy I]pi,(fc) — 
Yi^jj + Yi^jj{k) is independent of the frequency. By di- 
agonalizing the effective Hamiltonian /i^ -I- S(fc), one can 
explicitly carry out the frequency summation in Eq. (|13p . 
The parameters 61 and 62 are determined by 



'^1 = ^E["ll(^)-"44(fc)], (14) 

k 

^ ^E["22(fc)-"33(fc)]. (15) 
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FIG. 3: Equation for the particle-hole propagator D{k). The 
triangle denotes D{k). The effective interaction between par- 
ticles and holes is obtained by disconnecting a Green's func- 
tion line in the self-energy given in Fig. 2. 



So far, all the components of self-energy and parameters 
are self-consistently determined by Eqs. (|4|)- ([T5|) . The 
magnitude of 5i is larger than that of 62- To see it, 
consider temporarily the isolated bi and a2 atoms with- 
out Coulomb interaction. Since they are bonded by ti, 
their atomic degenerate states are split in two bonding- 
antibonding states with eigenvalue ztti. Therefore, the 
states of the bi and a2 sublattices contribute mostly to 
the eigenstates corresponding to the noninteracting en- 
ergy bands of overall energy separation ±ti from the zero 
energy. At low temperature, the lower band is occupied 
while the upper band is empty. On the other hand, the 
valence and conduction bands close to zero energy have 
eigenvectors which are composed predominantly of the 
linear combination of atomic states of the ai and b2 sub- 
lattices. The atoms of these two latter sublattices are 
the first to be affected by the Coulomb interaction, and 
they are subject to the most charge accumulation in the 
case of the GBISP. The two parameters Si and 62 are 
not independent, but are correlated through the Green's 
functions as described by Eqs. (g]) and (IT^ - p^ . We 
can chose 5i as the independent order parameter of the 
GBISP. 

To determine the GBISP phase boundary, that is the 
relation between the critical temperature Tq and the car- 
rier doping concentration S, we expand the self-energy 
and the Green's function to first order in the order pa- 
rameter Si. Let us define the matrix 



D{k) 



with 



d 

dSl 



[S(fc)-5^S*(fc)5]/2 



(16) 



^0 IN 

10 

10 

a oy 



(17) 



Notice that S*(fc) = S*(/c) (the transpose of E) since 
S^(fc) = S(fc). By this symmetry relation and by the 
definition in Eq. (fT6| . D{k) has the following structure: 



D 




D12 Di3 

D22 -Di3 

-£'22 -£'12 

-i^t3 -Dl^ -Diij 



(18) 



Therefore, only four elements Dn, D12, D13 and D22 
need to be determined. Under the mean-field approxi- 
mation, we have the following equation for D{k): 



k'XX' 



with 



J-fj.2 



dS2 
dSi 



(19) 

(20) 
(21) 



where G{k,iujeys are the normal-state Green's functions 
in which Si = S2 = 0. Again, the frequency summation in 
Eq. ((2T|) can be performed analytically. For the normal 
state, the Green's functions satisfy the relation G^i/ = 
Gpp with p, = 5- fj.. We therefore have /^^' = /l^^. By 
noting these relations, dS2/dSi can be expressed as 



(22) 



fcAA' 



Taking the partial derivative of Eq. ([14]) with respect to 
(5i, we obtain the condition for the phase transition. 



M 



^AY'(fc)i?AA'(fc) = l. 



(23) 



fcAA' 



Note that dn can be formally expressed as 



= ]g E Mil' (k) + "M2/2Y' MDxx' (k) 



fcAA' 



^ E V^Mf^^'{k)Dxyik), 



knXX' 



XX' 



where in the second line, the definition of u^i(2)i f^i) 



and D-^y{k) 



-D\i\{k) has been used. (The fac- 
tor 2 originates from the spin degeneracy.) Putting this 
result into Eq. (fT9| , one obtains the coupled linear equa- 
tions for ZJ's. The equations are diagrammatically shown 
in Fig. 3. The function D{k) is actually the particle-hole 
propagator. The effective interaction between particles 
and holes is the result of disconnecting a Green's func- 
tion in the self-energy as given in Fig. 2, by following the 
procedure explained in Fig. 3. Clearly, these equations 
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FIG. 4: Phase diagram of bilayer graphene in the SCHFA for 
coupling constant g — 1. The symbols are the numerical solu- 
tion for transition points. The dashed line is an extrapolation 
of the finite-temperature results to low temperature. 



are equivalent to solving the problem of a particle-hole 
propagator with a unity eigenvalue. 

Instead of solving the eigenvalue equations as given in 
Fig. 3, Z?(fc)'s can be determined more easily from Eqs. 
p9)) . ((20| and ([22]) by self-consistent iteration. For a 
given carrier doping concentration S, the transition tem- 
perature To can be found by gradually lowering tempera- 
ture T from a value higher than the critical one, and solv- 
ing the equations for Z?'s in the normal state at each step 
of the process. When the left-hand side of Eq. (1^ be- 
comes equal to 1, the critical temperature Tq is reached. 

To numerically solve Eqs. (g]), ^ and (Il9l)-(l23l), we 
need to carefully treat the convolution of the Coulomb 
interaction v^^^ (q) and the function n^^ (fc — q) as appear- 
ing in Eq. (|12p [and the similar one appearing in Eq. 
(1191) ] because v^i,{q) is singular aX q — Q. In Appendix 
A, we present an algorithm to deal with this problem. 

In Fig. 4, we show the phase diagram of the electron 
system in the 5 — T plane for coupling constant g = 1. 
At low temperature and low carrier doping, the system is 
in the GBISP. The transition temperature as a function 
of 5 is uniquely defined only at low 5 < 0.24 x 10"^. 
However, in the region 0.24 x lO""^ < 6 < 0.3 x lO"'^, 
each S corresponds to two transition temperatures. In 
the latter case, the phase boundary was determined by 
adjusting S for every given T. 

The numerical results for the order parameters di and 
62 as functions of T at (5 = for coupling constants 
g = 0.5, 1 and 1.8 are shown in Fig. 5. The SCHFA 
results are denoted as HF. We notice that \Si\ > \S2\, but 
62 is not negligibly small, which is different from what 
has been assumed in the two-band model. ^ From Fig. 5 
one can understand that the charge configuration at the 
four sites in the unit cell is (— |(52|, — 1<52|, |<^i|) [an- 



FIG. 5: Order parameters Si and ^2 as functions of temper- 
ature T at 5 = for coupling constants g = 0.5, 1 and 1.8. 
The symbols refer to numerical results. Circles and squares 
denoted as HF are the SCHFA results for 5i and ^2, respec- 
tively. The RRDA results denoted by triangles (Si) and in- 
verse triangles (52) are vanishingly small. 



other solution is — |(S2|, |<^2|, — The signs of Si 
and 62 are the opposite of each other because with such 
a charge distribution, the Coulomb interaction between 
sites a and b in the same plane is attractive and stabilizes 
the GBISP. It is also seen from Fig. 5 that the transi- 
tion temperature is higher for a system with stronger 
coupling. 

Our lattice model is different from both the two- and 
four-band continuum effective modelsi^'S'2£i2Z The two- 
and four-band continuum models are established under 
the consideration that the energy scale of quasiparti- 
cle spectral resonances involved in the problem is small 
with respect to a characteristic energy taken from bilayer 
graphene noninteracting band structure. For the four- 
band continuum model, the energy should be much less 
than the bandwidth of the tt orbitals of graphene. The 
two-band model for BLG is accurate only in the case in 
which the quasiparticle energy is much smaller than the 
gap ti. In the presence of long-range Coulomb interac- 
tion v{q), the energy transfer at small q is very large and 
the assumption for the validity of the two- and four-band 
continuum models becomes problematic. In this sense, 
the lattice model appears to be more reasonable. 

Another important difference between the lattice 
model and the two- and four-band continuum models 
relates to the valley physics in the Brillouin zone. In 
the two- and four-band continuum models, the two val- 
leys are independent of each other and the valley index 
is treated as an overall degeneracy index. On the con- 
trary, within the lattice model two states belonging to 
different valleys can be connected by nonzero intervalley 
Hamiltonian matrix elements. 
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FIG. 6: Additional part of the self-energy besides the Hartree- 
Fock terms. 



IV. SUPPRESSION OF THE GBISP IN THE 
RRDA 

The order parameters 5i and 82 so obtained by the 
SCHFA are overestimated because charge fluctuations 
have been ignored. We here reexamine the possibihty 
of the existence of the GBISP using the RRDA. 

Under the RRDA, besides the Hartree-Fock terms 
given in Fig. 2, the additional part of the self-energy, 
denoted by S'^(A:, iwf), is shown in Fig. 6. Each bubble 
in Fig. 6 is composed of two Green's functions G, rep- 
resenting the charge polarizability. In terms of G, the 
elements of ^'^{k, iloi) are expressed as 

T,''^^{k,iuJi) = -—^G^,^{k - q,iujt- ii'm)WI;^(q,iiym) 

where Vm is the bosonic Matsubara frequency, and 
W^^{q, ivrn) is an effective interaction. The matrix form 
of is given by 

W'{q, iv^) = [1 - v{q)x{q. WmT^v{q) - v{q) (24) 

with 

Xyivin^Wm) = —'^G^„{k,iuJi)G^f_,{k ~ q,iuji ~ iiy,n) 

k,f. 

and v[q) is the Fourier component (4x4 matrix) of the 
Coulomb interaction. The total self-energy is given by 

Y^^_,u{k, iujt) = Af_,S^iiy + S^^(fc) -I- S^^(fc, iuji). (25) 

The Green's function G is self-consistently determined 
and satisfies the microscopic conservation laws.^ 

Note that E'^ is a convolution of G and W^, and x is 
a convolution of two G's in momentum and frequency 
space. The easy way to calculate them is by Fourier 
transform. At low temperature, the summations index 
over the Matsubara frequencies should run up to a large 
frequency cutoff. To reduce the requirement of computer 
memory storage and accelerating the numerical compu- 
tation, the special algorithm of Ref. [25| can be used. 

The interaction W'^{q,Wrn) vanishes for m — > 00. For 
finite Vm , we need to carefully deal with the singularity at 
9 = 0. The Fourier transform W'^{q,iv,n) to W'^{r,ii>m) 
is discussed in Appendix B. 

At the ground state for T = 0, the Matsubara frequen- 
cies uji and Vm are treated as the continuous variables lo 



and V, respectively, and the summations over them are 
replaced by integrals, 



We computed the Green's function within the RRDA. 
The results for the order parameters 5i and 82 for g = 0.5, 
1 and 1.8 at (5 = are shown in Fig. 5 and compared 
with the SCHFA. The doping concentration chosen cor- 
responds to (5 = 0, for which the SCHFA transition tem- 
perature reaches its maximum. Though J = is the most 
favorable case for the GBISP predicted by the SCHFA, 
the two order parameters are substantially suppressed 
by CDF in the RRDA; the magnitude of the two param- 
eters is three orders smaller than that of the SCHFA. 
For inspecting the GBISP ordering at low temperature, 
the RRDA calculation for g = 1.8 is performed down to 
T = 0. From the numerical results, we conclude that 
there is no the GBISP in systems with 5 < 1.8 under the 
RRDA. 

We also computed the Green's function within RPA, 
in which the polarizability xilj'^^m) in W'^{q,ivm) [see 
Eq. is replaced by the polarizability for noninter- 

acting electrons. At (5 = 0, similar to the RRDA, the 
parameters 5i and 82 so obtained are vanishingly small. 
In the RRDA, a replacement of the bare Coulomb inter- 
action in the Hartree term by the screened one is prohib- 
ited because where the ring diagrams are equivalent to a 
self-energy insertion to the Green's function to be renor- 
malized. When performing the RPA calculation, we also 
need to keep the bare Coulomb interaction in the Hartree 
term. 

The reason for the suppression of the GBISP un- 
der the RRDA is that the exchange interaction is sig- 
nificantly weakened by the screening due to electronic 
charge-density fluctuations while the Hartree term op- 
posing the charge transfer between the two layers^ is not 
changed. At low temperature, in a wide range of Matsub- 
ara frequencies, the exchange interaction is short-ranged 
and weakened and does not favor the GBISP transition. 
We point out that the suppression of the GBISP here is 
not due to prohibition by the Mermin- Wagner theorem<^ 
The theorem applies to a system with continuous sym- 
metry; if the symmetry were broken, there would be 
a logarithmically diverging number of long-wavelength 
collective fluctuations accompanying the excitations on 
top of the broken symmetry ground state of the two- 
dimensional system. In the present case, the inversion is 
a discrete symmetry operation, and there is no diverging 
long-wavelength collective fluctuation arising from the 
breaking of such a symmetry. 
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V. SUMMARY 

In summary, we have studied the physics of interact- 
ing electrons in bilayer graphene using the lattice model. 
The possibility of the existence of a GBISP at low tem- 
perature and low-carrier-doping concentration is reinves- 
tigated with both the SCHFA and the RRDA. The lat- 
ter approach takes into account the charge density fluc- 
tuations beyond the SCHFA or the mean-field approxi- 
mation. Under the RRDA, the exchange interaction is 
weakened substantially, and the existence of a GBISP 
becomes unsustainable. We have also presented the nu- 
merical method for dealing with convolution of a singular 
Coulomb interaction and the Green's function on the lat- 
tice model. This numerical method should be usable for 
solving problems in many-particle systems. 
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Appendix A: Calculation of the momentum- 
space convolution of Coulomb interaction v^rith a 
smooth function for a lattice model 

For solving problems of two-dimensional electron sys- 
tem in the presence of long-range Coulomb interaction, 
we sometimes need to deal with the convolution 



Cik) 



q) 



(26) 



where the q summation runs over the first Brillouin zone, 
V{q) is the Coulomb interaction, and /(fc) is a smooth 
function of k. On a lattice, an analytical expression for 
V{q) is not available but its long- wavelength behavior is 
known. For the honeycomb lattice, it is given by Eq. 
(fTO|) . V{q) can be divided into long-range and short- 
range parts. For the honeycomb lattice under considera- 
tion, define 



\Qn + Q\ 



exp(-ao|Q„ + Q|) (27) 



where c — 27re^/S'oe is the same factor that appeared in 
Eq. pop . Qn is the reciprocal lattice vector, Q = Mq 
is as given in the text, and Gq is an auxiliary parameter. 
By taking ao = 2a, the summation in Eq. (|27l) converges 
quickly and only a few terms need to be summed up. 
Clearly, «'((?) represents a long-range interaction. With 
'y'(g), V{q) can be written as 



where (q) is so defined by the equation and is the short- 
range part of V{q). Note that both w'(q) and w'*(g) are 
periodic functions of q. Equation (|26|) now is given by 



The first integral in Eq. ([291) can be safely performed by 
Fourier transform. In the second integral, the singularity 
appears at q = 0. To find out an auxiliary function for 
this integral, we pay attention to the expanding form of 

f{k-q) 



f{k -q) ^ f{k) - qxfx{k) - qyfyik) 



(30) 



where fx(y) (k) ~ df{k) / dk^(^yj . Define two auxiliary func- 
tions Vx{q) and Vy{q) by 

Vx{y){q) = 2^ -T — exp(-ao|g„ + Ql), 

where v^i^y^(q) is periodic and odd under q — > — g. The 
second integral in Eq. (|29l) can be written as 

±J2^^(q)f{k~q) = l-J2{v\k^k')[fik')-f{k)] 



M 

fc' 

+Vxik--k')Mk) 

+ Vy{k-k')fy{k)} 

+f{k)vHr)\r=o- 



(31) 



Now, there is no singularity in the integrand in Eq. 
(ini). The leading term of v\k - k')[f{k') - f{k)] as 
fc' — fc is proportional to the derivative of / multiplied 
with a sign factor since v''{k — k') oc l/\M{k' — k)\. 
This leading term varies discontinuously at fc' = k. 
The discontinuity is canceled out by the remaining term 
Vx{k — k')fx{k) + Vy{k — k')fy{k). As a result, the in- 
tegrand is a smooth function. The integral can then be 
carried out numerically. The last term stems from the 
introduction of the auxiliary functions to the integrand. 
The value u'(r)|r=o is given by 



(32) 



which can be calculated explicitly. Replace g-summation 
by 



1 

M 



dQ 



(33) 



V(q)=v\q)+v\q) 



(28) 



where Sq — ^/ia? /2 is the area of the unit cell of the hon- 
eycomb lattice as appearing in the text, and BZ means 
the integral is performed over the first Brillouin zone. 
The combination of the integration over BZ and the Q^- 
summation in the definition of v^{q) equals the integra- 
tion of the function cexp(— aoQ)/Q over the total space 
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of Q, 



dQ 



(27r)2 Q 



exp(-ao(3) 



e^/ao. 



(34) 



The function /(fc) is assumed to be smooth here. How- 
ever, for calculating the Fock exchange self-energy, f{k) 
corresponds to the distribution function and varies dra- 
matically at the Fermi surface at low temperature. In 
this case, extremely dense grids in a momentum regime 
covering the Fermi surface should be used to denote the 
variation of f{k). 

The term Vx{k — k')fx{k) + Vy{k — k')fy{k) was in- 
troduced in the right-hand side of Eq. (I5T]) in order to 
smooth the integrand. Because Vx{q) and Vy{q) are pe- 
riodic and odd functions of q, the contribution from the 
integral of Vx{k — k')fx{k) -\- Vy{k — k')fy{k) is zero. At 
T = 0, there is a discontinuity in f(k) at the Fermi sur- 
face and its derivatives fx{k) and fy{k) are 6 functions. 
Therefore, the use of this term at T = is unworthy. At 
T = 0, this term should be removed, keeping the discon- 
tinuity in the integrand. The cost is to use dense grids 
near the Fermi surface to ensure the accuracy of the re- 
sult. 

Appendix B: Fourier transform of the screening 
potential W^{q,ivm) 

To take the Fourier transform of W'^{q^ivm) given by 
Eq. (j24l) from momentum space to real space, we first 
pay attention to its singularity at q = 0. For small 
because x{q,ii^m) is finite, the singularity exists only in 
the second term v{q) on the right-hand side of Eq. (IMl) . 
Its real space form is known as that given by Eq. ^ for 
its elements. However, at large z/^, because xC?,*^'™) is 
vanishingly small, there is also a singularity in the first 
term on right-hand side of Eq. (1^^ and both terms can- 
cel with each other. We need a systematic numerical 
scheme for the Fourier transform at any z/„j. 

Note that in the limit q 0, v{q) — > VQ{q)A with 
vo{q) = c/Q (again c = 27re^/S'oe and Q = \Mq\) and 



and WmiQ) so defined by Eq. (p5)) . By observing this 
asymptotic form, we take the auxiliary function for the 
Fourier transform as 



where ao again is a parameter for fast convergence of 
the summation over the reciprocal lattice vectors Qn- 
The Fourier transform of W^qjivm) is separated into 
two parts, [W''{q,Wrn) - Wa{q)A] and Wa{q)A. There is 
no singularity in the first one and it can be safely trans- 
formed by numerical computation. For the second one, 
Wa{q) is transformed as 



Wa{r) 



dq 



So 
So 



BZ (2^) 
dQ 



BZ (27r 
dQ 



^Wa{q) exp{iq- r) 
^Wa{q) exp{iQ-R) 



(27r)^ 



Jo 



2tt 



dQ'-^^^MQR) (38) 
Q + am 



where the first line is the definition with q and r given 
in the quadrilateral coordinate system, the second line 
converts q to Q = Mq and R = (M*)"V (with M* 
the transpose of M) in the orthogonal systems with 
dq — dQ/\M\ — ^/3dQ/2, the third line comes from the 
definition of Wa{q) given by Eq. (I57| . the last line is 
obtained after the azimuthal integration, and Jo{QR) is 
the Bessel function. Now the singularity in the integrand 
exists only when = 0, but also appears in the 
front factor and the integral vanishes. However, for large 
-R, JoiQR) oscillates rapidly with Q. By observing the 
large-Qi? behavior of Jo{QR), we choose the auxiliary 
function^^ 




In the same limit, we have 



Q{Q + a„ 

W,niQ)A, 



-A 



with 



Ja{z) = 



1 



TTZ + 1 
[1- 



{[1 



8(nz + ly 
I cos(z)} 



(39) 



and separate Jo{QR) to Jo{QR) — Ja{QR) and Ja{QR)- 
(-35) By replacing MQR) with Jo{QR)- JaIqR) in Eq. ([38]), 
the integral can be accurately carried out by simple nu- 
merical method. The remaining integral with Jo{QR) 
replaced by Ja{QR) can be performed using Filon's 
(36) method. 
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